home *** CD-ROM | disk | FTP | other *** search
/ SGI Hot Mix 17 / Hot Mix 17.iso / HM17_SGI / research / lib / comfit.pro < prev    next >
Text File  |  1997-07-08  |  6KB  |  157 lines
  1. ;$Id: comfit.pro,v 1.7 1997/01/15 03:11:50 ali Exp $
  2. ;
  3. ; Copyright (c) 1994-1997, Research Systems, Inc.  All rights reserved.
  4. ;       Unauthorized reproduction prohibited.
  5. ;+
  6. ; NAME:
  7. ;       COMFIT
  8. ;
  9. ; PURPOSE:
  10. ;       This function fits the paired data {X(i), Y(i)} to one of six common
  11. ;       types of appoximating models using a gradient-expansion least-squares
  12. ;       method. The result is a vector containing the model parameters.
  13. ;
  14. ; CATEGORY:
  15. ;       Statistics.
  16. ;
  17. ; CALLING SEQUENCE:
  18. ;       Result = COMFIT(X, Y, A)
  19. ;
  20. ; INPUTS:
  21. ;       X:    An n-element vector of type integer, float or double.
  22. ;
  23. ;       Y:    An n-element vector of type integer, float or double.
  24. ;
  25. ;       A:    A vector of initial estimates for each model parameter. The length
  26. ;             of this vector depends upon the type of model selected.
  27. ;
  28. ; KEYWORD PARAMETERS:
  29. ;       EXPONENTIAL:    If set to a non-zero value, the parameters of the 
  30. ;                       exponential model are computed. Y = a0  * a1^x + a2.
  31. ;
  32. ;         GEOMETRIC:    If set to a non-zero value, the parameters of the
  33. ;                       geometric model are computed.  Y = a0 * x^a1 + a2.
  34. ;    
  35. ;          GOMPERTZ:    If set to a non-zero value, the parameters of the
  36. ;                       Gompertz model are computed.  Y = a0 * a1^(a2*x) + a3.
  37. ;
  38. ;        HYPERBOLIC:    If set to a non-zero value, the parameters of the
  39. ;                       hyperbolic model are computed.  Y = 1./(a0 + a1*x)
  40. ;
  41. ;          LOGISTIC:    If set to a non-zero value, the parameters of the
  42. ;                       logistic model are computed.  Y = 1./(a0 * a1^x + a2)
  43. ;
  44. ;         LOGSQUARE:    If set to a non-zero value, the parameters of the
  45. ;                       logsquare model are computed.
  46. ;            Y = a0 + a1*alog10(x) + a2 * alog10(x)^2
  47. ;
  48. ;             SIGMA:    Use this keyword to specify a named variable which 
  49. ;                       returns a vector of standard deviations for the computed
  50. ;                       model parameters.
  51. ;                
  52. ;           WEIGHTS:    An n-element vector of weights. If the WEIGHTS vector 
  53. ;                       is not specified, an n-element vector of 1.0s is used.
  54. ;
  55. ;              YFIT:    Use this keyword to specify a named variable which 
  56. ;                       returns n-element vector of y-data corresponding to the 
  57. ;                       computed model parameters.
  58. ;
  59. ; EXAMPLE:
  60. ;       Define two n-element vectors of paired data.
  61. ;         x = [2.27, 15.01, 34.74, 36.01, 43.65, 50.02, 53.84, 58.30, 62.12, $
  62. ;             64.66, 71.66, 79.94, 85.67, 114.95]
  63. ;         y = [5.16, 22.63, 34.36, 34.92, 37.98, 40.22, 41.46, 42.81, 43.91, $
  64. ;             44.62, 46.44, 48.43, 49.70, 55.31]
  65. ;       Define a 3-element vector of initial estimates for the logsquare model.
  66. ;         a = [1.5, 1.5, 1.5]
  67. ;       Compute the model parameters of the logsquare model, a(0), a(1), & a(2).
  68. ;         result = comfit(x, y, a, sigma = sigma, yfit = yfit, /logsquare)
  69. ;       The result should be the 3-element vector:
  70. ;         [1.42494, 7.21900, 9.18794]
  71. ;
  72. ; REFERENCE:
  73. ;       APPLIED STATISTICS (third edition)
  74. ;       J. Neter, W. Wasserman, G.A. Whitmore
  75. ;       ISBN 0-205-10328-6
  76. ;
  77. ; MODIFICATION HISTORY:
  78. ;       Written by:  GGS, RSI, September 1994
  79. ;-
  80.  
  81. pro exp_func, x, a, f, pder
  82.   f = a[0] * a[1]^x + a[2]
  83.   if n_params() ge 4 then $
  84.     pder = [[a[1]^x], [a[0] * x * a[1]^(x-1)], [replicate(1., n_elements(x))]]
  85. end
  86.  
  87. pro geo_func, x, a, f, pder
  88.   f = a[0] * x^a[1] + a[2]
  89.   if n_params() ge 4 then $
  90.     pder = [[x^a[1]], [a[0] * alog(x) * x^a[1]], [replicate(1., n_elements(x))]]
  91. end
  92.  
  93. pro gom_func, x, a, f, pder
  94.   f = a[0] * a[1]^(a[2]*x) + a[3]
  95.   if n_params() ge 4 then $
  96.     pder = [[a[1]^(a[2]*x)], [a[0] * a[2] * x * a[1]^(a[2]*x-1)], $
  97.     [a[0] * x * alog(a[1]) * a[1]^(a[2]*x)], [replicate(1., n_elements(x))]]
  98. end
  99.  
  100. pro hyp_func, x, a, f, pder
  101.   f = 1.0 / (a[0] + a[1] * x)
  102.   if n_params() ge 4 then $
  103.     pder = [[-1.0 / (a[0] + a[1] * x)^2], [-x / (a[0] + a[1] * x)^2]]
  104. end
  105.  
  106. pro log_func, x, a, f, pder
  107.   f = 1.0 / (a[0] * a[1]^x + a[2])
  108.   if n_params() ge 4 then begin
  109.     denom =  -1.0/(a[0] * a[1]^x + a[2])^2
  110.     pder = [[a[1]^x*denom], [a[0] * x * a[1]^(x-1)*denom], [denom]]
  111.     endif
  112. end
  113.  
  114. pro logsq_func, x, a, f, pder
  115.   b = alog10(x)
  116.   b2 = b^2
  117.   f = a[0] + a[1] * b + a[2] * b2
  118.   if n_params() ge 4 then $
  119.     pder = [[replicate(1., n_elements(x))], [b], [b2]]
  120. end
  121.  
  122. function comfit, x, y, a, weights = weights, sigma = sigma, yfit = yfit, $
  123.                        exponential = exponential, geometric = geometric, $
  124.                        gompertz = gompertz, hyperbolic = hyperbolic, $
  125.                        logistic = logistic, logsquare = logsquare
  126.  
  127.   on_error, 2
  128.  
  129.   fcn_names = ['exp_func', 'geo_func', 'gom_func', 'hyp_func', 'log_func', $
  130.         'logsq_func']
  131.   fcn_npar = [ 3, 3, 4, 2, 3, 3]
  132.  
  133.   nx = n_elements(x)
  134.   wx = n_elements(weights)
  135.   if nx ne n_elements(y) then $
  136.     message, 'x and y must be vectors of equal length.'
  137.  
  138.   if wx eq 0 then weights = replicate(1.0, nx) $
  139.   else if wx ne nx then $
  140.     message, 'x and weights must be vectors of equal length.'
  141.  
  142.   a1 = A            ;Copy initial guess
  143.   if keyword_set(exponential) then i = 0 $
  144.   else if keyword_set(geometric) then i = 1 $
  145.   else if keyword_set(gompertz) then i = 2 $
  146.   else if keyword_set(hyperbolic) then i = 3 $
  147.   else if keyword_set(logistic) then i = 4 $
  148.   else if keyword_set(logsquare) then i = 5 $
  149.   else message, 'Type of model must be supplied as a keyword parameter.'
  150.  
  151.   if n_elements(a1) ne fcn_npar[i] then $
  152.     message, 'A must be supplied as a '+strtrim(fcn_npar[i], 2) + $
  153.          '-element initial guess vector.'
  154.   yfit = curvefit(x, y, weights, a1, sigma, function_name = fcn_names[i])
  155.   return, a1
  156. end
  157.